<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of demhmc3</title>
  <meta name="keywords" content="demhmc3">
  <meta name="description" content="DEMHMC3 Demonstrate Bayesian regression with Hybrid Monte Carlo sampling.">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html &copy; 2003 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="../../m2html.css">
</head>
<body>
<a name="_top"></a>
<div><a href="../../menu.html">Home</a> &gt;  <a href="#">ReBEL-0.2.7</a> &gt; <a href="#">netlab</a> &gt; demhmc3.m</div>

<!--<table width="100%"><tr><td align="left"><a href="../../menu.html"><img alt="<" border="0" src="../../left.png">&nbsp;Master index</a></td>
<td align="right"><a href="menu.html">Index for .\ReBEL-0.2.7\netlab&nbsp;<img alt=">" border="0" src="../../right.png"></a></td></tr></table>-->

<h1>demhmc3
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="box"><strong>DEMHMC3 Demonstrate Bayesian regression with Hybrid Monte Carlo sampling.</strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="box"><strong>This is a script file. </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="fragment"><pre class="comment">DEMHMC3 Demonstrate Bayesian regression with Hybrid Monte Carlo sampling.

    Description
    The problem consists of one input variable X and one target variable
    T with data generated by sampling X at equal intervals and then
    generating target data by computing SIN(2*PI*X) and adding Gaussian
    noise. The model is a 2-layer network with linear outputs, and the
    hybrid Monte Carlo algorithm (with persistence) is used to sample
    from the posterior distribution of the weights.  The graph shows the
    underlying function, 300 samples from the function given by the
    posterior distribution of the weights, and the average prediction
    (weighted by the posterior probabilities).

    See also
    <a href="demhmc2.html" class="code" title="">DEMHMC2</a>, <a href="hmc.html" class="code" title="function [samples, energies, diagn] = hmc(f, x, options, gradf, varargin)">HMC</a>, <a href="mlp.html" class="code" title="function net = mlp(nin, nhidden, nout, outfunc, prior, beta)">MLP</a>, <a href="mlperr.html" class="code" title="function [e, edata, eprior, mse] = mlperr(net, x, t)">MLPERR</a>, <a href="mlpgrad.html" class="code" title="function [g, gdata, gprior] = mlpgrad(net, x, t)">MLPGRAD</a></pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../../matlabicon.gif)">
<li><a href="hmc.html" class="code" title="function [samples, energies, diagn] = hmc(f, x, options, gradf, varargin)">hmc</a>	HMC	Hybrid Monte Carlo sampling.</li><li><a href="mlp.html" class="code" title="function net = mlp(nin, nhidden, nout, outfunc, prior, beta)">mlp</a>	MLP	Create a 2-layer feedforward network.</li><li><a href="mlpfwd.html" class="code" title="function [y, z, a] = mlpfwd(net, x)">mlpfwd</a>	MLPFWD	Forward propagation through 2-layer network.</li><li><a href="mlpinit.html" class="code" title="function net = mlpinit(net, prior)">mlpinit</a>	MLPINIT Initialise the weights in a 2-layer feedforward network.</li><li><a href="mlppak.html" class="code" title="function w = mlppak(net)">mlppak</a>	MLPPAK	Combines weights and biases into one weights vector.</li><li><a href="mlpunpak.html" class="code" title="function net = mlpunpak(net, w)">mlpunpak</a>	MLPUNPAK Separates weights vector into weight and bias matrices.</li></ul>
This function is called by:
<ul style="list-style-image:url(../../matlabicon.gif)">
<li><a href="demnlab.html" class="code" title="function demnlab(action);">demnlab</a>	DEMNLAB A front-end Graphical User Interface to the demos</li></ul>
<!-- crossreference -->


<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="fragment"><pre>0001 <span class="comment">%DEMHMC3 Demonstrate Bayesian regression with Hybrid Monte Carlo sampling.</span>
0002 <span class="comment">%</span>
0003 <span class="comment">%    Description</span>
0004 <span class="comment">%    The problem consists of one input variable X and one target variable</span>
0005 <span class="comment">%    T with data generated by sampling X at equal intervals and then</span>
0006 <span class="comment">%    generating target data by computing SIN(2*PI*X) and adding Gaussian</span>
0007 <span class="comment">%    noise. The model is a 2-layer network with linear outputs, and the</span>
0008 <span class="comment">%    hybrid Monte Carlo algorithm (with persistence) is used to sample</span>
0009 <span class="comment">%    from the posterior distribution of the weights.  The graph shows the</span>
0010 <span class="comment">%    underlying function, 300 samples from the function given by the</span>
0011 <span class="comment">%    posterior distribution of the weights, and the average prediction</span>
0012 <span class="comment">%    (weighted by the posterior probabilities).</span>
0013 <span class="comment">%</span>
0014 <span class="comment">%    See also</span>
0015 <span class="comment">%    DEMHMC2, HMC, MLP, MLPERR, MLPGRAD</span>
0016 <span class="comment">%</span>
0017 
0018 <span class="comment">%    Copyright (c) Ian T Nabney (1996-2001)</span>
0019 
0020 
0021 <span class="comment">% Generate the matrix of inputs x and targets t.</span>
0022 ndata = 20;                     <span class="comment">% Number of data points.</span>
0023 noise = 0.1;                    <span class="comment">% Standard deviation of noise distribution.</span>
0024 nin = 1;                        <span class="comment">% Number of inputs.</span>
0025 nout = 1;                       <span class="comment">% Number of outputs.</span>
0026 
0027 seed = 42;                    <span class="comment">% Seed for random number generators.</span>
0028 randn(<span class="string">'state'</span>, seed);
0029 rand(<span class="string">'state'</span>, seed);
0030 
0031 x = 0.25 + 0.1*randn(ndata, nin);
0032 t = sin(2*pi*x) + noise*randn(size(x));
0033 
0034 clc
0035 disp(<span class="string">'This demonstration illustrates the use of the hybrid Monte Carlo'</span>)
0036 disp(<span class="string">'algorithm to sample from the posterior weight distribution of a'</span>)
0037 disp(<span class="string">'multi-layer perceptron.'</span>)
0038 disp(<span class="string">' '</span>)
0039 disp(<span class="string">'A regression problem is used, with the one-dimensional data drawn'</span>)
0040 disp(<span class="string">'from a noisy sine function.  The x values are sampled from a normal'</span>)
0041 disp(<span class="string">'distribution with mean 0.25 and variance 0.01.'</span>)
0042 disp(<span class="string">' '</span>)
0043 disp(<span class="string">'First we initialise the network.'</span>)
0044 disp(<span class="string">' '</span>)
0045 disp(<span class="string">'Press any key to continue.'</span>)
0046 pause
0047 
0048 <span class="comment">% Set up network parameters.</span>
0049 nhidden = 5;            <span class="comment">% Number of hidden units.</span>
0050 alpha = 0.001;                  <span class="comment">% Coefficient of weight-decay prior.</span>
0051 beta = 100.0;            <span class="comment">% Coefficient of data error.</span>
0052 
0053 <span class="comment">% Create and initialize network model.</span>
0054 
0055 <span class="comment">% Initialise weights reasonably close to 0</span>
0056 net = <a href="mlp.html" class="code" title="function net = mlp(nin, nhidden, nout, outfunc, prior, beta)">mlp</a>(nin, nhidden, nout, <span class="string">'linear'</span>, alpha, beta);
0057 net = <a href="mlpinit.html" class="code" title="function net = mlpinit(net, prior)">mlpinit</a>(net, 10);
0058 
0059 clc
0060 disp(<span class="string">'Next we take 100 samples from the posterior distribution.  The first'</span>)
0061 disp(<span class="string">'300 samples at the start of the chain are omitted.  As persistence'</span>)
0062 disp(<span class="string">'is used, the momentum has a small random component added at each step.'</span>)
0063 disp(<span class="string">'10 iterations are used at each step (compared with 100 in demhmc2).'</span>)
0064 disp(<span class="string">'The step size is 0.005 (compared with 0.002).'</span>)
0065 disp(<span class="string">'The new state is accepted if the threshold'</span>)
0066 disp(<span class="string">'value is greater than a random number between 0 and 1.'</span>)
0067 disp(<span class="string">' '</span>)
0068 disp(<span class="string">'Negative step numbers indicate samples discarded from the start of the'</span>)
0069 disp(<span class="string">'chain.'</span>)
0070 disp(<span class="string">' '</span>)
0071 disp(<span class="string">'Press any key to continue.'</span>)
0072 pause
0073 
0074 <span class="comment">% Set up vector of options for hybrid Monte Carlo.</span>
0075 nsamples = 100;        <span class="comment">% Number of retained samples.</span>
0076 
0077 options = foptions;     <span class="comment">% Default options vector.</span>
0078 options(1) = 1;        <span class="comment">% Switch on diagnostics.</span>
0079 options(5) = 1;        <span class="comment">% Use persistence</span>
0080 options(7) = 10;    <span class="comment">% Number of steps in trajectory.</span>
0081 options(14) = nsamples;    <span class="comment">% Number of Monte Carlo samples returned.</span>
0082 options(15) = 300;    <span class="comment">% Number of samples omitted at start of chain.</span>
0083 options(17) = 0.95;    <span class="comment">% Alpha value in persistence</span>
0084 options(18) = 0.005;    <span class="comment">% Step size.</span>
0085 
0086 w = <a href="mlppak.html" class="code" title="function w = mlppak(net)">mlppak</a>(net);
0087 <span class="comment">% Initialise HMC</span>
0088 <a href="hmc.html" class="code" title="function [samples, energies, diagn] = hmc(f, x, options, gradf, varargin)">hmc</a>(<span class="string">'state'</span>, 42);
0089 [samples, energies] = <a href="hmc.html" class="code" title="function [samples, energies, diagn] = hmc(f, x, options, gradf, varargin)">hmc</a>(<span class="string">'neterr'</span>, w, options, <span class="string">'netgrad'</span>, net, x, t);
0090 
0091 clc
0092 disp(<span class="string">'The plot shows the underlying noise free function, the 100 samples'</span>)
0093 disp(<span class="string">'produced from the MLP, and their average as a Monte Carlo estimate'</span>)
0094 disp(<span class="string">'of the true posterior average.'</span>)
0095 disp(<span class="string">' '</span>)
0096 disp(<span class="string">'Press any key to continue.'</span>)
0097 pause
0098 
0099 nplot = 300;
0100 plotvals = [0 : 1/(nplot - 1) : 1]';
0101 pred = zeros(size(plotvals));
0102 fh1 = figure;
0103 hold on
0104 <span class="keyword">for</span> k = 1:nsamples
0105   w2 = samples(k,:);
0106   net2 = <a href="mlpunpak.html" class="code" title="function net = mlpunpak(net, w)">mlpunpak</a>(net, w2);
0107   y = <a href="mlpfwd.html" class="code" title="function [y, z, a] = mlpfwd(net, x)">mlpfwd</a>(net2, plotvals);
0108   <span class="comment">% Sum predictions</span>
0109   pred = pred + y;
0110   h4 = plot(plotvals, y, <span class="string">'-r'</span>, <span class="string">'LineWidth'</span>, 1);
0111 <span class="keyword">end</span>
0112 pred = pred./nsamples;
0113 <span class="comment">% Plot data</span>
0114 h1 = plot(x, t, <span class="string">'ob'</span>, <span class="string">'LineWidth'</span>, 2, <span class="string">'MarkerFaceColor'</span>, <span class="string">'blue'</span>);
0115 axis([0 1 -3 3])
0116 
0117 <span class="comment">% Plot function</span>
0118 [fx, fy] = fplot(<span class="string">'sin(2*pi*x)'</span>, [0 1], <span class="string">'--g'</span>);
0119 h2 = plot(fx, fy, <span class="string">'--g'</span>, <span class="string">'LineWidth'</span>, 2);
0120 set(gca, <span class="string">'box'</span>, <span class="string">'on'</span>);
0121 
0122 <span class="comment">% Plot averaged prediction</span>
0123 h3 = plot(plotvals, pred, <span class="string">'-c'</span>, <span class="string">'LineWidth'</span>, 2);
0124 
0125 lstrings = char(<span class="string">'Data'</span>, <span class="string">'Function'</span>, <span class="string">'Prediction'</span>, <span class="string">'Samples'</span>);
0126 legend([h1 h2 h3 h4], lstrings, 3);
0127 hold off
0128 
0129 disp(<span class="string">'Note how the predictions become much further from the true function'</span>)
0130 disp(<span class="string">'away from the region of high data density.'</span>)
0131 disp(<span class="string">' '</span>)
0132 disp(<span class="string">'Press any key to exit.'</span>)
0133 pause
0134 close(fh1);
0135 clear all;</pre></div>
<hr><address>Generated on Tue 26-Sep-2006 10:36:21 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/">m2html</a></strong> &copy; 2003</address>
</body>
</html>